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ABSTRACT 

Models of accretion discs and their associated outflows often incorporate as- 
sumptions of axisymmetry and symmetry across the disc plane. However, for 
turbulent discs these symmetries only apply to averaged quantities and do not 
apply locally. The local asymmetries can induce local imbalances in outflow 
power across the disc mid-plane, which can in turn induce local tilting torques. 
Here we calculate the effect of the resulting stochastic torques on disc annuli that 
are a consequence of standard mean field accretion disc models. The torques in- 
duce a random walk of the vector perpendicular to the plane of each averaged 
annulus. This random walk is characterized by a radially dependent diffusion 
coefficient which we calculate for small angle tilt. We use the coefficient to cal- 
culate a radially dependent time scale for annular tilt and associated jet wobble. 
The wobble time depends on the square of the wander angle so the age of a given 
system determines the maximum wobble angle. We apply this to examples of 
blazars, young stellar objects and binary engines of pre-planetary nebulae and 
microquasars. It is noteworthy that for an averaging time t w ~ 3 days, we esti- 
mate a wobble angle for jets in SS433 of 9 ~ 0.8 degrees, not inconsistent with 
observational data. In general the non-periodic nature of the stochastic wobble 
could distinguish it from faster periodic jet precession. 



Introduction 



Winds and jets are commonly associated with accr e tion disc engines, and standard 
models of discs (IShakura fc Sunvaevl Il973l; iPringla Il98ll ; iFrank et al.l 120021 ) and outflows 



(jBlandford fc Payne 



1982 



Pelletier fc Pudritall992l ) typically impose axisymmetry and re- 



flection symmetry across the disc plane. However, theory and observation suggest that for 
a wide range of sou rces, angu l ar mo mentum transport is likely mediated at least in part 



by turbulence (e.g. iTerquem I (120021 )). This immediately implies that the assumptions of 
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axisymmetry apply at most only to mean flow quantities, averaged over sufficient time or 
space (e.g. Balbus et al. 1994). The local violation of symmetries for a turbulent disc, leads 
to an intrinsic variability (Blackman 1998). 

We suggest that an additional consequence of the fluctuations may be a local "rocket 
effect." When a wind is expelled from the surface of a disc it exerts a force on that surface 
which compensates for the momentum of the outflow. In idealized models of symmetric 
accretion discs the forces exerted on the disc by the winds on either side of the disc cancel 
out. In reality, although these forces must cancel globally if the disc center of mass is 
stationary, turbulence induced asymmetry makes it unlikely that these forces cancel out 
locally. Thus, different portions of the disc can be displaced above or below the initial 
mean symmetry plane. Here we investigate how disc turbulence can cause a local imbalance 
between the winds emanating from the top and bottom by averaging the local effect over 
azimuth. The net effect is to induce a stochastic wander of the vector perpendicular to the 
plane determined by the azimuthal average of the asymmetry for each eddy at a given radius. 
We model this as Brownian motion, and calculate the associated diffusion coefficient. 

In Sec. 2 we discuss the needed basics of Brownian motion. In Sec. 3 we then derive 
the equation for the aforementioned stochastic diffusion coefficient for tilt as a function of 
disc parameters. In Sec. 4 we apply this diffusion coefficient to the wobbling of a thin 
accretion disc and a stochastic wander of the associated jet axis, assumed to be parallel to 
the disc's primary angular velocity vector. We apply the formalism to four astrophysical 
engine classes (blazars, young stellar objects (YSOs), planetary nebulae, microquasars) and 
estimate the wander angles for which the wobble could be observable in these systems. As 
a particular example, we compare the crudely predicted wobble to that observed in SS433 
and find agreement to order of magnitude for a plausible choice of disc parameters. The 
stochastic nature of the tilt distinguishes it from that of an ordered precession We conclude 
in Sec. 5. 



2. Review of Brownian Motion 



Consider the random walk on a two-dimensional square lattice of points (I Chandr asekhar 



19431 ). We label a point on the lattice by coordinates (6, 0), and the spacing between adjacent 
lattice points by A9. We also break time into discrete intervals At. We use these coordinates 
in anticipation of our application to the random walk of a unit vector tip on the surface of 
a sphere, locally approximated as a plane for small angles. 



At time t — At, a point particle located at (9, (f>) moves by time t to one of the 4 points 
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(9 + A8, <p),{9- A9, 0), (9, + A9), {9, - AO) with equal probability. If the particle moves 
into the lattice point (9, 0) at time t, it would have come from one of the four adjacent 
lattice points in the previous time step. Therefore, if P(6, 0, t) is the probability the particle 
is located at (9, 0) at time t, we have 

P(6, 0, t) = \{P{6 + A9, 0, t-At) + P(9 - A9, 0, t - At) 

+P{9, + A9,t - At) + P(6, - A9,t - At)) ^ ' 

This allows us to write a difference equation for P(9,(p,t). Using ([1]) we have 

P(9,(p,t)-P(9,(j),t-At) 
At 

where 



D(A(9,<P) + B(9,<f>)) (2) 



P(6 + A9, 0, t - At) - 2P(9, 0, t-At) + P{8 - A9, 0, t - At) 
A(9, 0) = — , (3) 



P(g, + Afl, t - At) - 2P(6, 0, t - At) + P(9, - Ag, t - At) 

0) = m2 , (4) 

and 

4At v ; 

Equation (j2j) is a discrete version of the two-dimensional diffusion equation that characterizes 
a random walk with diffusion coefficient D. Although we will use the discrete version in 
this paper, a continuum diffusion equation results by taking the continuum limit for which 
A# -> and At -> 0. In this limit, A -> and B -> 0. If we then define D = 
li m At,A6»-^o(^^) 2 /^? we have: 

D h^ + ir^) (6) 



dt V. 

which is a two-dimensional Cartesian diffusion equation in the local angle coordinates. 



3. Derivation of the Wobble Diffusion Coefficient 

To derive D for a disc, we first consider a small volume element at radius r with az- 
imuthal and radial scales equal to that of a turbulent eddy and with vertical scale equal to 
the disc thickness. We define dL as the net mechanical wind luminosity emanating from this 
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volume element, and v as the associated mean outflow speed at the chosen radius. We posit 
that the winds emanating from either side of thin annulus are not exactly in momentum bal- 
ance at all times due to turbulent fluctuations. We denote wind mechanical luminosities for 
a given eddy emanating from the top and bottom of the disc as dL top and dL bot respectively, 
such that dL = dL top + dLb ot . Let e top and e bo t denote unit vectors in the direction of wind 
flow on the top and bottom. The net force on the volume element from the two winds is 

7-i-i ,tti dL top dLbot /dL 

at — at t0 p + at bot — e top -\ e bot — e — , { ( ) 

V V V 

where e' = {e top dL t o V + e bot dLhot) / dL, and characterizes the direction of wind imbalance. 
Since we have assumed that the origin of the force imbalance (turbulence) is intrinsically 
random, e' will be a random vector dependent on both time and position. Fluctuations 
around both the mean outflow speed and mass outflow rate can be incorporated into e' given 
its definition; the explicitly presence of the mean velocity in (JTj) does not preclude this. 

The force on a disc volume element results from a "rocket effect" created by non- 
vanishing last term in (JTj). To calculate the associated torque, we assume that e' is a poloidal 
vector (e' ■ <p = 0) for every volume element. The associated torque then lies in the initial 
disc plane and is given by 

dr = (r x dF) = (r x e')— = (f x e')-dL . (8) 

v v 

Summing the torques for all volume elements at radius r leaves a net torque on an annulus 
of radial thickness dr at radius r given by 

v-^ , r r dL 

dT net = } (r x e )-dL = e — — dr , (9) 
^— ' v v dr 

dV&E{r) 

where we have defined E(r) = {all volume elements at radius r}, and e = ^(r x e'). Note 
that since e' is poloidal, r x e' lies in the plane of the disc, and therefore so does e, being a 
sum of such vectors. Since dr lies in the plane of the disc, so does dT net . By summing over all 
the volume elements in the annulus we have removed the azimuthal dependence, providing 
what is needed to characterize the response of an axisymmetric mean field annulus. This is 
an annulus of a "standard" axisymmetric radially dependent accretion disc (which itself is 
axisymmetric precisely because of the azimuthal average over turbulent fluctuations). 

We can now calculate the rate of change of angular momentum of this annulus associated 
with tilting from its initial disc symmetry plane. This tilt introduces an angular momentum 
vector in the initial disc plane that bisects the annulus, and is approximately perpendicular 
to the underlying Keplerian rotation vector. The component of the moment of inertia tensor 
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appropriate for the tilt axis is dl = ^r 2 dm(r), where the mass of the annulus is dm(r) = 
fj,(r)dr and /i(r) is the mass per unit radius. The rate of change in angular momentum that 
characterizes the tilt is then 

duj r 2 duj r 2 duj 

d3 = ^ dI = J^ dm= ^H^ r)dr > (10) 

where uj is the orbital angular momentum vector that would reflect a fixed-plane orbital 
motion of disc material in the absence of rocket torques. Using Newton's second law dJ = 
dT net . It follows that 

-^(r) = e r -^. (11) 
2 dt v dr 

If we assume that e varies on a characteristic time scale At, then ffTTl) implies 

2 dL 

Aoj = — — eAt. 12 

rvfiyr) dr 

Eq. (fl2j) measures the tilt of the plane the annulus. Since dT ne t is in the plane of the 
annulus, the torque has the effect of rotating uj but not changing its magnitude. If the 
angular deflection is small in a time At, then AO ~ where AO is the change in the angle 
of uj during the time At. We can then appeal to the Cartesian lattice caluclation of Sec. 2 
to obtain the small angle diffusion coefficient: Combining the expression for AO just derived 
with ffl2l and ([5]) we obtain 



r 2 ^(rfv 2 uo 2 \dr J r 2 /j,(rfv 2 uj 2 \ dr J y/N 

where e e d = N l l 2 e represents the characteristic contribution from each of the N eddies at a 
given radius to e, and r e ^ = r/N 1 ^ 2 is the characteristic time scale for evolution of a given 
eddy. Here N = 2nr /l e d is the number of eddies with scale l e d at radius r. 

This diffusion coefficient characterizes the random walk of the orbital velocity vector of 
the mean field annulus due to tilt from the stochastic rocket force. Since the magnitude of 
oj never changes in our approximation, the tilt implies that the orbital angular vector moves 
on the surface of a sphere as the annulus wobbles. For small tilt angles, this surface can be 
locally approximated by a plane and we are justified in using the Cartesian formulation of 
Sec. 2, despite the global differences between spheres and planes. 

The derivation and use of Eqn. (fT5|) implies that each radius has its own wobble dif- 
fusion coefficient. This does not violate angular momentum conservation because annuli 
are (at least) viscously coupled. If all of the material that enters the disc at large radii is 
coplanar, then conservation of total angular momentum implies that the total integrated 
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angular momentum about any axis perpendicular to that of the initial ui must vanish at all 
times. This means that if some inner mean annulus of gas tilts in one direction, there will 
be a corresponding counter-tilt elsewhere. Turbulent viscosity allows the radial distribution 
of tilts to be non-trivial and because jet power is likely dominated by inner annuli, the tilt of 
the inner-most annuli are particularly relevant for predicting jet wobble because they dom- 
inate the accretion power. In what follows, we assume that the axis of any annular section 
of a jet tracks the orbital axis of the corresponding disc annulus to which that jet section is 
anchored. Thus the jet wobble directly tracks the disc wobble. 

4. Application to Astrophysical Engines 

Here we apply the above formulae to thin accretion discs and jet wobble. 



4.1. Determination of Stochastic Wobble times for Discs and Jets 

For a thin, non-self-gravitating disc in hydrostatic equilibrium 

uo = uo K = v'GM/r 3 ~ c s /H ~ tj (14) 
where c s is the sound speed, M is the central ob j ect m ass, and H is the density scale 



height. The Shakura-Sunyaev (jShakura &: Sunyaevl (119731 )) viscosity v in hydrostatic equi- 



librium then satisfies v = ac s H ~ v^ d /t e d ~ v^ d H/c s , where a is the viscosity parameter. 
Then (e.g. Blackman 1998) v e d ~ c^^Cg and 

l ed = a 1/2 H, (15) 



so 



27r r . , 

N = ^W (16 > 



The disc accretion rate can be modeled as M a (r) = M (r/r ) s (e.g. iBlandford fc Begelman 



( 119991 )). where M is the accretion rate at the outer edge of the disc, r Q is the radius of the 
outer edge, and < s < 1 is a parameter. The total mass outflow rate from r Q to r is then 
M{r) = M Q [1 — (r/r ) s }. Taking the derivative, we obtain the outflow mass loss rate by a 
small annulus of width dr at radius r to be 

dM , „> s ( r V" 1 , 

-dr = -M — — dr. (17) 



dr r Q \r a 
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The luminosity of the gas ejected by the wind in the annulus is dL 
outflow speed satisfies 

v ~ v esc = y/2GM/r, 
so in combination with (IT7j) we have 



dL • s ( r 
— = M o7 — - 



^v 2 \dM\, where the 



dr 



2r Q \r 



(18) 



(19) 



Using dH]), (USD, (USD, (HID and f[TOT) in (JI3) gives 



DM 



-i 



-1 



Af 



10 18 g/s 



rad 2 /s, (20) 



where r g = GM/c 2 . For a standard Shakura-Sunyaev disc model (IShakura fc Sunyaevill973l ) 



supplemented by our radially dependent accretion rate M = M {r/r ) s we have 



/j(r) = 5.1xlO n a- 4/5 



7/10 



10 18 g/s 



7s/10 



M 

Mm 



1/2 



1/4 



1- ^ 



r \ 1/2 



1 7/20 



g/cm, 
(21) 



and 



H 



9 x 10- 3 a^ 10 



M 



-1/4 



3/20 



10 18 g/s 



3s/20 



9/8 



r \ 1/2' 



3/20 



Using ( 1211) and ( T22l) in ( 1201) . the time required to wander an angle 9 radians is 



CM 



4D 



1.6 x 10 4 a~ 9 / 5 (-2- 



/ M 
V 10 18 g/e 



-27 
40 



17 27a 
16 40 



M 



17/8 / \ 27s/40 
ro 



(22) 

i-m 1/2 

(23) 



5/8 



sec. 



For stochastic wobble to be observable t w (8) must at least be less than the age of a 
system r age . A further constraint follows from the standard accretion paradigm in which 
accretion supplies the energy for a turbulent viscosity which, together with the outflow, 
transport angular momentum and sustain the accretion. As discussed at the end of Sec. 3, 
the turbulent viscosity is also a means of transporting the angular momentum associated 
with tilt between annuli such that the tilt of inner annuli is compensated by the opposing 
tilt of outer annuli. For stochastic wobble to produce a net mean tilt on an annulus at an 
inner radius, the viscous time scale there T u (r) should be short compared to the tilt time, so 
that the compensating angular momentum opposite to that associated with the tilt can be 
transported outward. Note that even though material also accretes on a viscous time, the 
resupply of material at a given radius would also maintain the stochastic equilibrium on time 
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scales short compared to the secular wobble time. Provided there is a supply of mass, we 
therefore consider viscosity to maintain the conditions for the effect to accrue, rather than 
washing out the effect at a given radius. 

Summarizing the above conditions to be checked, we have 

r v (r) < t w (r) < r age , (24) 
where, from (T2"2"j) . the viscous accretion time scale for r > 3r g is given by 




and 

r age = Max r v (r )/-^- (26) 

where M c is the companion mass. The second value in brackets is used when the accretion 
disc is being supplied by a companion, while the first (the viscous time at the outer disc 
radius) is used when there is no companion-fed resupply. 

For later use in evaluating ( 1241) . we have from ( 123)) and ( 125)) 



m ~ ,7 x ^ (£)' (*)-*"* (27) 

9 17 27s trr/o 27s 5 



Note that the left side of (f28|) employs r for the numerator and r Q in the denominator on 
the left hand side as that ratio represents a comparison between the wobble time at a given 
radius with the accretion time at the outer radius. Eq. (127)) compares the two time scales at 
the same radius. In the next section we estimate (123)) and check condition (EMI) for example 
classes of accretion engines. 



4.2. Application to Astrophysical Jets 

To find the characteristic jet/disc wobble times for various systems, we now estimate 
( 123)) for active galactic nuclei (AGN) young stellar objects (YSO), microquasars, planetary 
nebulae (PNe) and active galactic nuclei (AGN). Since the wander time scale depends on the 
square of the wander angle, we can determine the wander angle for which the observability 
condition (T25)) is satisfied. Presently, the quantities s and e e d, and a are not tightly pinned 
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by either theory or simulation so we maintain their dependence explicitly. Fiducial values 
will be taken to be e ec i = 0.5, and a mass-loss power law index s = 0.2, so their product can 
be scaled to 0.01. We also scale the results with a. 

Consider first AGN blazars (e.g. 3C273; Paltani and Turler 2005). We use M ~ 
1O 9 M and Mq ~ lOM /yr. Then taking tq = 10 G r g and r = r g ~ 10 14 cm, we find that 
t w ~ 3.1 x 10 11 (f§f)~ 2 (fi) ~ 9/5 ( 003 e rad ) 2 yr. For these parameters the condition 9 < 0.03 
corresponds to satisfying the second inequality in (1241 . namely that Eq. (1281) is < 1. However, 
observable duty cycles of AGN jets seem to last only ~ 10 8 yr, which would requires a wobble 
angle 9 < 100" . For this latter range of wobble angles, both inequalities in ff24|) are satisfied, 
using the first value on the right of ff26l) . 

For YSOs (Hartmann 1998), we take M = 1M , M ~ lO~ 5 M /yr, r Q = 100 AU. At 
r = 7 x 10 7 r 9 = 10 13 cm, Eq. ([23]) then gives t w ~ 4.6 (fgf)" 2 (j§)" 9/5 (o^) 2 yr. Using 
(|2T|) . the first inequality of ([M]) is easily satisfied. For the second inequality of (I2"4"[) we use 
the first value on the right of (G2|. From (EHD we then find ~ 0.3 (^) 2 satisfied 

for wander angles < O.Olrad. 

Accretors with binary companions supply accretion for times much longer than the 
viscous time scale at the outer radius. This makes the second equality in (T21]) easier to satisfy 
for large wobble angles in those systems. Two examples of such are accretors are pre-PNe 
engines and microquasars. For pre-PNe engines, an accretion disc can form around the post- 
ascending giant branch (AGB) white dwarf core as the result of common envelope evolution 
(Nordhaus & Blackman 2006). We take M = O.6M , M ~ lO" 5 M /yr, r = 10 n cm. 
Here the engine accretor is a proto- white dwarf of radius ~ 10 9 cm or r ~ 10 4 r 9 . At 
r = 10% ~ 10 10 cm, Eq. QZ® then gives t w ~ 5.5 (fgf)" 2 (^)" 9/5 (5^3) 2 yr. Again 
for this case the first inequality in (j24p is satisfied from (T27j): At the inner disc radius 
t w / T p ~ 700 ( 03 6 ' rad ) • In this case however, the second equality in tf2fi$ is more easily 
satisfied as we can use the second value in ( |26|) : The pre-PNe outflows last T age ~ 10 3 yr and 
the same accretion disc, supplied by a very low mass star, or very high mass planet, may in 
fact be continuously supplying outflow power even during the 10 4 yr PNe phase (Blackman 
et al. 2001). Both inequalities of (|24|) can then be satisfied for 9 < 0.8 rad. 

For microquasar examples we consider Cygnus X-l (Young et al. 2001) and SS433 (Katz 
& Piran 1982; Begelman et al. 2006; Blundell et al. 2005;2007). For the former, we take 
M = 1OM , M ~ 10 19 M o /yr, and r = W00r g . At r = r g ~ 1.5 x 10 6 cm, Eq. ([23]) then 

gives t w ~ 0.03 (fffi) 2 ((^) 9 ^ 5 (oaTTid) 2 V 1 - ^ or tms case ^ ne ^ rst inequality in fT24l) 
is easily satisfied from fT271) . The second inequality in ff24|) allows use of the second value in 
fT26|) . The companion mass to e.g. Cyg X-l is at least 2OM so for accretion at lO~ 7 M /yr, 
accretion can be powered long after the wobble time. 
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For SS433 (e.g. Begelman et al. 2006) M = 10M Q , and we take the inner radius 
r ~ r g ~ 1.5 x 10 6 cm. The accretion rate that makes it to the inner region, relevant for the 
jets, is M in ~ 2 x 10 19 g/ s. This is to be distinguished from the super-Eddington accretion rate 
of M out > 10 3 M in at the disc circularization radius of r circ ~ 10 12 cm. Begelman et al. (2006) 
argue that the much of the mass is taken out by a wind at r ~ ^gj 2 ~ 10 5 r 9 ~ 10 11 cm. We 
therefore take r Q ~ 10 5 r 9 for the outer radius of the disc supplying mass to the inner jet and 
use M in for the relevant accretion rate. We then obtain t w ~ 0.04 (^gf) 2 (^) 9 ^ 5 ( 03 6> rad ) 2 
yr. For SS433, like Cygnus X-l, the first inequality in (|2~4|) is easily satisfied from (|2"T|) and 
a companion mass of ~ 2OM (based on crudely averaging recent mass estimates of Lopez 
et al. 2006 and Hillwig and Gies 2008) allows use of the second value in fT26l) for the second 
inequality in ff24^) . This implies that accretion will be fueled long after the wobble time. 
Again, we have not considered relativistic effects. 

We can compare the stochastic wander just estimated for SS433 with the data of Blundell 
et al. (2005,2007), who monitored fluctuations in the speed and direction of the inner jets. 
Table 1 of Blundell et al. (2007) shows that when averaged over several day time scales, 
wander angles of ~ 0.5 — 0.7 degrees were obtained. This is roughly consistent with our 

ft \ 1/2 

result above: Converting units, out result for SS433 can be written 9 ~ 1.7 ( 15d a ys ) deg. 

There is likely a maximum time scale t w ^ max which is the value of t w above which the 
stochastic wander is tempered by some negative feedback from forces not considered in the 
present work. Thus t w < t w ^ max if the wander is to be calculated with our approximations. 
One limit on t w comes from the small wander angle approximation that allowed our use of 
a 2-D Cartesian approximation. However, for a precessing and nodding system like SS433 
(e.g. Fabrika 2004), the 162 day precession time scale or the 6.3 day nodding time scale 
may provide the maxium t w appropriate for computing untempered stochastic wander. We 
do not determine t w ^ max here or study large angle wobble but note that a persistent feature 
of stochastic wander is its non-periodicity, which distinguishes it from systematic precession 
induced by binary companions (e.g. Terquem et al. 1999). 



5. Conclusions 

Turbulence in accretion discs violates local axisymmetry and reflection symmetry even 
if these properties hold in some average mean sense. As such, this symmetry violation 
likely also applies to outflows emanating from the disc. Due to a rocket-like effect, any 
outflow power asymmetry locally displaces the disc material in the direction opposite to the 
net outflow. By azimuthally averaging outflows with locally reflection asymmetric powers, 
we have calculated the radially dependent small-angle diffusion coefficient of the angular 
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momentum vector for disc annuli. The diffusion of this vector represents a stochastic tilting 
of the annulus defined by the averaging, and implies a disc warping and jet wobble. For the 
latter, the most relevant tilt is that associated inner disc radii. We have calculated estimated 
the stochastic wobble angle vs. time scale for jets and applied the relation to a handful of 
astrophysical engine classes to exemplify crude application of the paradigm. The numbers 
are not inconsistent with measured values in SS433. We have not considered the nonlinear 
saturation of the stochastic wobble in the present work. 

EGB acknowledges support from NSF grant AST-0406799, NASA grant ATP04-0000- 
0016, and NORDITA during the Turbulence and Dynamos workshop 2008. 
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